Solving the general 3-D safety factor by combining Sarma’s idea with the assumption of normal stress distribution over the slip surface

This study proposes a method for determining 3-D limit equilibrium solutions. The method, inspired by Sarma, introduces the horizontal seismic coefficient as a slope failure parameter and implements a modification of the normal stress over the slip surface. Four equilibrium equations are used to solve the problem without compromising the accuracy of the calculations: three force equilibrium equations in the x, y, and z directions and a moment equilibrium equation in the vertical (z) direction. The reliable factor of safety can be determined by calculating the minimum value of the horizontal seismic coefficient. Furthermore, we analyzed several typical examples of symmetric and asymmetric slopes, finding good consistency with the existing literature. This consistency indicates the reliability of the factor of safety we obtained. The proposed method is favored due to its straightforward principle, convenient operation, fast convergence, and ease of programming.


Introduction
Slope stability is a classical problem in soil mechanics, and it is ubiquitous in the field of geotechnical engineering. Existing methods include the limit analysis method [1,2], the limit equilibrium method [3], the strength reduction method [4], and the variational analysis [5]. Most of which exist because there is no established unifying method for slope stability. Two-dimensional (2-D) slope stability analysis methods are most commonly used among engineers due to their simplicity. These 2-D methods are based on simplified assumptions, and in this way, they reduce the dimension of the problem from 3-D to 2-D. For a known slip surface, 2-D methods are based on the limit equilibrium theory that usually employs the method of slices, including Bishop [6], Morgenstern and Price [7], Spencer [8], Janbu [9], Sarma [10], Fredlund and Krahn [11], and so on.
The complexity of slip surfaces in 3D has rendered 2D methods incapable of yielding precise results. Accordingly, tentative development of 3D methods has taken place, including those by Xing [12], Hungr et al. [13], and Huang and Tsai [14]. For a given slip surface, the factor of safety is predominantly calculated under partial equilibrium conditions using non-rigorous 3D methods, such as those developed by Lam and Fredlund [15], Huang et al. [16], and a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 Chen and Zhu [17]. Efforts have been made by Hovland [18], Hungr et al. [13], and Boutrup and Lovell [19] to extend 2D models, originally developed by Fellenius [3], Bishop [6], and Janbu [9], to their 3D counterparts. However, such extensions are associated with assumptions about inter-column forces, which may lead to divergent results in iterative calculations.
The safety factor can be obtained if the distribution of the normal stress can be determined in the slope stability analysis. The assumed distribution of the normal stress acting on the slip surface is proposed by Bell [20], which is a function with two parameters to be determined. Zhu and Lee [21] assume that the distribution of the normal stress acting on the slip surface is composed of an initial function and a modified function, the latter of which was assumed to be a linear interpolation function with two parameters to be determined. This method does not need to assume the inter-column forces, thereby compensating for flaws in the integration of calculations in conventional methods. Zhu et al. [22] extend the assumption of the normal stress to the 3D limit equilibrium, only considering the vertical force equilibrium and the moment equilibrium by rotating axis for symmetric slopes. Moreover, they establish a 3-D calculation method that satisfies different equilibrium conditions for the slope stability based on this concept [23][24][25]. Zhu and Qian [26] suggest that the modified normal stress distribution along the sliding surface should be non-negative. The methods proposed by Zhu and Qian [26] and Lu and Zhu [27] need to solve the strict limit equilibrium solution from a nonlinear system that contains six variables.
Our team presents a simplified method for calculating the factor of safety of a slope with a three-dimensional (3-D) slip surface. This method satisfies all six equilibrium conditions for the 3-D sliding body. By incorporating Sarma's concept, we introduce a clever trick that employs the horizontal seismic force coefficient to linearize the equations [28]. Through this, we obtain a rigorous limit equilibrium solution for 3-D slope stability. However, there are cases where negative normal stress may occur over the surface edge, necessitating modifications. This process is lengthy and complicated, it is difficult for engineers to control and apply. Therefore, we propose utilizing four equilibrium equations to address the problem without compromising the accuracy of the calculation. It involves disregarding two minor equilibrium conditions, the moments about the x-and z-axis, rendering it easily applicable to engineering calculations. The efficacy of the proposed method has been evaluated by applying it to a selection of both symmetric and asymmetric classic three-dimensional slope examples. The results have demonstrated that the method is both quick and reliable. Fig 1(A) proposed a 3-D collapsed slope, it is assumed that the direction of the slip surface is opposite to the x-axis, and the sliding mass is discretized into n columns. The slope and slip surface are described by functions s(x,y) and g(x,y), respectively. There are three forces acting on a typical discretized column, including: (i) the weight of the column, w(x,y), which numerically equal gðx; yÞ½sðx; yÞ À gðx; yÞ�, where γ(x,y) is the unit soil weight, τ(x,y) is the shear stress at the base, (ii) the normal stress at the base, σ(x,y), and (iii) the horizontal seismic coefficient, K c , as shown in Fig 1(B). The discretization widths of the column in the x-and y-directions are defined as dx and dy, respectively. α x and α y are the horizontal inclination angles of the slip surface along the x-and y-directions, respectively. According to the geometric of the column information, we have:

Basic assumptions
cosg z ¼ 1 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi D ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi s 02 where γ z is the angle between the normal of the slip surface and the z-axis. s 0 x and s 0 y are the tangent values that the angle between the column base and x-and y-axis directions. In the threedimensional slope stability analysis, it is necessary to analyze and calculate the directions of both normal stress (σ) and shear stress (τ). The cosine direction of the normal stress over the slip surface can be represented by ðn x s ; n y s ; n z s Þ, where: The direction of the potential sliding mass is generally parallel to the x-axis, while the direction of shear stress is perpendicular to the direction of normal stress. Therefore, the cosines of the shear stress τ can be written as ðn x t ; n y t ; n z t Þ, where: Where D ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi 1 þ s 0 x 2 p . As shown in Fig 2, there are many rectangles when the slip surface of a small column is projected to the x-y plane. Each rectangle has a base area of dA, which can be expressed as: The initial normal stress over the slip surface is obtained by applying Hovland method [18]: where w is the weight of the column. In this paper, the proposed normal stress distribution over the slip surface is described by Eq (13), while the modification function is shown in Eq (14), both of which are from Zhu and Qian [26]: where σ(x,y) is the total normal stress distribution when the sliding mass is in the limit equilibrium state. Based on the limit equilibrium methods, the proposed method considers the sliding mass as a rigid mass in the limit equilibrium state, and it satisfies four equations of limit equilibrium according to the assumption of normal stress. Herein, the modification function is defined as where λ 1 , λ 2 and λ 3 are unknown parameters.

The limit condition
To ensure the consistency with the limit equilibrium approach and the validity of the Mohr-Coulomb criterion, the shear strength, τ, is defined as: where c 0 and ψ are the soil cohesion and the friction coefficient, respectively; ψ = tanφ 0 , where φ is the effective angle of internal friction; c u = c 0 −uψ, where u is the pore water pressure and F s is the slope safety factor.

Equilibrium equations of the failure mass
In a three-dimensional sliding body, the proposed approach satisfies the force equilibrium in all three directions and the moment equilibrium in the vertical direction: Substituting Eq (16) into Eqs (17A)-(17D) respectively yields: ∬s n x s s À n z s x Eqs (18A)-(18D) can be abbreviated as:

Numerical solution for the factor of safety
K c is set as a known physical parameter to discuss its effect on the slope stability analysis. The proposed method determines the factor of safety, F s for a slope in a different way than conventional methods, because it considers two to-be-known parameters, namely K c and F s . Usually, F s is assumed to be 1, while K c required to produce this F s is solved as an unknown. The linear equations of the four variables can be obtained by substituting Eq (13) into Eq (19), yielding: Eq (20) can be rewritten in the matrix form: When i = 1, 2, 3 and j = 1, 2, 3, we have: When i = 1,2,3,4 and j = 4, we have: The matrix can quickly solve the four equations, and their formulas are listed in Table 1.
As is indicated in Eq (20), K c exists in a linear equation system and can, through the matrix, solve Eq (22C).
To calculate F s for a given slip surface, K c is needed to implement several iterations by changing the value of F s . Eq (23) is used to search for the minimum K c and the corresponding Table 1. Calculation parameters for three-dimensional slope stability analysis.

PLOS ONE
Combining Sarma's idea for 3-D safety factor F s :

Solution for the safety factor
This study aims to find the minimum K c that can determine F s . Specifically, F s is negatively correlated with K c ; that is, a larger horizontal seismic force is needed to cause the slope sliding for a more stable slope [29,30]. Therefore, K c would be very small if the slope is in a very stable state. There are several steps to compute F s using the present framework, as shown in Fig 3: 1. Divide the failure mass into n columns.
2. Calculate the assumed initial normal stress over the slip surface σ 0 and its modification ξ i .

Calculate the values of
4. Assume an initial iteration value of the safety factor F s as 1, which is expressed as F 1 .
5. Compute the value of K 1 and the second value of safety factor F 2 .
6. Replace F 1 by F 2 and return to calculate K 2 .
7. Calculate the value of F 3 according to Eq (26).
8. Use the value of F 3 to calculate the value of K 3 .
9. Calculate the difference ΔF 1 between F 3 and F 1 , and the difference ΔF 2 , between F 2 and F 1 .
10. Replace F 2 by F 3 and K 3 by K 2 when ΔF 1 is smaller than ΔF 2 ; otherwise, replace F 2 by F 1 , K 2 by K 1 , F 3 by F 1 , and K 3 by K 2 .
The iteration procedure is to repeatedly change the values of F s to get different value of K c until |K 3 −K 2 |<ε that indicates K 3 to be close to zero (Part 2 in Fig 3). Hence, the final value of K c corresponds to the actual factor of safety of the three-dimensional slope. Due to the linear relationship between K c and λ, only several times of iterations are needed to calculate the actual value of F s .

Numerical examples and verification
A computer program is developed in MATLAB to illustrate the 3D shape of the failure surface and calculate its safety factor. Four examples are performed to evaluate the feasibility of this program for stability analysis of slopes with ellipsoid and asymmetric slip surfaces.

Example 1: A ellipsoid slip surface with weak layer
The proposed method is first validated in 2 cases using the 3D example dataset provided by Xing [12]. Case 1 has a failure surface that is assumed to be a symmetric ellipsoid, while Case 2 has a failure surface that is composed of an ellipsoid surface and a weak layer. The 3D slope profile used in Example 1 is presented in Fig 4. They have been studied by many researchers, including Xing [12], Hungr et al. [13], Huang and Tsai [14], Chen and Zhu [17], Zhu and Qian [26], Lu and Zhu [27], Zhu et al. [28] and Chen et al. [31]. The cross-section and parameters are shown in Fig 5, including a slope height (H) As shown in Table 2, the proposed method obtained similar values of F s on the composite failure surface with those by previous methods for Case 1. The safety factor calculated by the proposed method is 2.142, compared to the 2.122 value given by Xing [12]. Both methods satisfy force equilibrium in three directions and moment equilibrium in the vertical direction, signifying that the set of equilibrium conditions is similar for both. However, the result from

PLOS ONE
Combining Sarma's idea for 3-D safety factor our proposed method is larger than that of Xing [12], which can be attributed to the omission of the intercolumn force during our calculation. The difference between the method used in this study and the Rigorous 3-D methods used by Zhu and Qian [26] and Lu and Zhu [27]

PLOS ONE
Combining Sarma's idea for 3-D safety factor is.0.75% and 0.05%, respectively, which is essentially negligible. This fact validates the feasibility of the proposed method in solving for the safety factor using four equilibrium equations. The high similarity between the results of the proposed method and Zhu and Qian [28] may stem from the same assumptions used in these two methods. Compared with other methods such as those by Hungr et al. [13], Huang and Tsai [14], Chen and Zhu [17], and Chen et al. [31], the errors are 1.17%, 3.41%, 3.17%, and 2.1%, respectively. These errors fall within an acceptable range, demonstrating that the results obtained by this method are in line with those from previous methods and thus are acceptable. In Case 2, as shown in Table 2, the safety factor calculated by the proposed method is 1.662. An 7.76% difference exists between the results of the proposed method and those from Xing [12]. This discrepancy could be attributed to the shear force on the two side faces of all columns along the slip direction, which was ignored in our proposed method. Similarly, the results from our method are consistent with those of Lu and Zhu [27]. The discrepancies between our method and the methods of Zhu and Qian [26] and Zhu et al. [28] are 0.12% and 0.06%, respectively. The high similarity in safety factors can be attributed to the assumption employed regarding the normal stress distribution along the slip surface. The difference with Hungr et al. [13], Huang and Tsai [14], Chen and Zhu [17], and Chen et al. [31] may be caused by the following factors: (1) Different number of the discrete columns; (2) Different assumption of inter-column force; (3) Different assumptions of normal stress distribution.
As shown in Fig 6, the safety factor of the slope increases as the horizontal seismic coefficient decreases in all iteration steps.
Various columns were used to discretize the failure mass to examine the effects of the number of columns on solution accuracy. Fig 7 illustrates the grid sensitivity. The safety factor ranges from 2.139 to 2.215 in Case 1, while it varies between 1.654 and 1.675 in Case 2. The safety factor gradually stabilizes as the number of columns increases. Therefore, it can be indicated that the safety factor converged to a certain value as the columns increased.
Verification of the modified normal stress. Figs 8 and 9 show the distributions of initial normal stress and the modified normal stress over the slip surface of Case 1 and Case 2, respectively. The modified normal stresses over the slip surface are positive, smooth, and continuous. Therefore, the results are reasonable.

Example 2: A homogenous slope from Lu and Zhu [27]
Example 2 is based on the dataset provided by Lu and Zhu [27]. Fig 10 schematically illustrates that a homogenous slope with a height of H = 40m, and strength parameters of soil c 0 = 30 kPa, φ 0 = 30 kPa, and γ = 22 kN/m 3 . The slip surface is an ellipsoid described by the equation: where B is the length of the slip surface in x direction (Fig 11). For L = 1,2,3,6 and 8 times the height of the slope, the proposed method's safety factor F s is calculated by the summaries in Table 3. This example intends to discuss the effect of the length of slip surface on the safety factor and compare the effect to that of different three-dimensional methods. The trend of the safety factor becomes more stable as the L/H ratio increases, and it shows good agreement with all methods from Table 3. Fig 12 indicates that the distribution of the modified normal stress over the slip surface is positive. By confirming that the value of the modified normal stress is non-negative, these results verify that the current method is reliable and acceptable.

PLOS ONE
Combining Sarma's idea for 3-D safety factor

Example 3: A circular slip surface adapted from Xing [12]
The slope height H = 12m and the slope ratio is 1:2; There are two cases in this example, (1) Case 3 with cohesion of c 0 = 20kPa, internal angle of φ 0 = 20˚, and unit weight γ = 18.8 kN/m 3 , and (ii) Case 4 with c 0 = 0, and φ 0 = 20˚. The cross-section and parameters of the spherical slope are shown in Fig 13. The slip surface can be numerically described as: ðx À 15:1Þ 2 þ y 2 þ ðz À 19:1Þ 2 ¼ 24 2 ð26Þ Table 4 presents a comparison of calculated Fs values using different methods. The discrepancies between the current method and that of Liu [32] are about 3.57% and 1.47% for Case 3 and 4, respectively. This discrepancy is due to Liu [32] employing a rigorous limit equilibrium method, which yields relatively conservative results. The proposed method, when compared to the method of Zhu et al. [28], shows only minor differences of 0.08% and 0.06% for Cases 3 and 4, respectively. This demonstrates the proposed method's effectiveness in analyzing slope

PLOS ONE
Combining Sarma's idea for 3-D safety factor stability for the circular slip surface, as indicated by the minimal error rates. As shown in Table 5, K c decreases as F s increases, in both Cases 3 and 4. Reliable F s is determined when K c approaches to zero.

Example 4: Wedge failure adapted by Hoek and Bray [33]
The stability of wedge in rock mechanics is a typical issue concerning three-dimensional limit equilibrium. The calculated examples involve both symmetric and asymmetric geometric wedge shapes. To test the validity of the proposed method, the results obtained by the proposed method are compared with those by Zhu et al. [28], Guo et al. [34], Jiang and Zhou [35], and Su and Shao [36], respectively. The geometric shape diagram for symmetric and  As to the symmetric wedge, c 0 = 0.02MPa and φ 0 = 20˚, while as to the asymmetric wedge, c 0 = 0.05MPa, and φ 0 = 30˚. The rock unit weight is 26 kN/m 3 . The parameter and geometry information using in the method are listed in Table 6. As shown in Table 7, there exist minor (and practically negligible) differences between the values calculated using the proposed method and those calculated by other methods. These

PLOS ONE
Combining Sarma's idea for 3-D safety factor discrepancies could be due to the imperfect replication of slip surface geometry. For an asymmetric wedge, Su and Shao [36] present significantly different results, as they take into account the effect of Poisson's ratio.
A strong correlation exists between K c and F s , as shown in Fig 18. F s = 1 is used as the initial safety factor for six steps of iteration. This case used 1 to start iteration of the initial safety factor, and ended the program through six steps. The values of instantaneous K c decrease with  the increase of F s values at step 2 and keep decreasing until the safety factor is stable. This phenomenon is a result of the fact that the horizontal seismic coefficient has a great influence on the safety factor.   Fig 19 displays a comparison between the distribution of the initial normal stress and the distribution of the modified normal stress for asymmetry and symmetry wedge failure. The modified normal stress sketches show that the modified normal stress is non-negative, verifying the rationality of the results.

Discussion of the advantages
This paper eliminated two moment balance equations considered relatively minor in engineering to simplify the problem-solving process. Only three force equilibrium equations in the (x, y, z) directions and a one-moment equilibrium equation in the (z) direction were considered. This method streamlines the calculation process without compromising computational accuracy. Table 7. A comparison of the safety factors for the wedge solution.

Methods
The factor of safety, F s

Symmetric Asymmetric
Zhu et al. [ Compared with the rigorous 3-D method, this method simplified the formulas of computation and avoided the negative normal stress at the edges of the slip surface. We found that disregarding the two secondary moment equilibrium conditions does not affect the accuracy of the results. It points out that the error of Fs is within the acceptable level. As shown in Table 8, compared to the rigorous method, the error with Case 1 is 0.04%, with Case 2 is 0.06%, with